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^ ■ We show in this paper by using the Wang-Landau flat-histogram Monte Carlo 

^ , method that the phase transition in the XY stacked triangular antiferromagnet is 

clearly of first-order, confirming results from latest Monte Carlo simulation and 

'^ ■ from a nonperturbative renormalization group, putting an end to a long-standing 

PACS numbers: 75.10.-b General theory and models of magnetic ordering ; 75.40.Mg Nu- 
merical simulation studies 
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O ' Effects of the frustration in spin systems have been extensively investigated during the 

last 30 years. Frustrated spin systems are shown to have unusual properties such as large 

(^ I ground state (GS) degeneracy, additional GS symmetries, successive phase transitions with 

complicated nature. Frustrated systems still challenge theoretical and experimental meth- 

%> ! ods. For recent reviews, the reader is referred to Ref. [l|. 



Let us confine our discussion on the nature of the phase transition in strongly frustrated 
spin systems. Since the nature of the phase transition depends on the symmetry and the 
dimension of the system, we have to examine first its GS properties. Of course, the nature 
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of the order parameter defined according to the system symmetry determines the properties 
of the phase transition. One of the most studied systems is the stacked triangular anti- 
ferromagnet (STA): the antiferromagnetic (AF) interaction between nearest-neighbor (NN) 
spins on the triangular lattice causes a very strong frustration. It is impossible^ to fully 
satisfy simultaneously the three AF bond interactions on each equilateral triangle. The GS 
configuration of both Heisenberg and XY models is the well-known 120-degree structure. 

The phase transition in strongly frustrated spin systems is a subject of intensive inves- 
tigations in the last 20 years. Theoretically, these systems are excellent testing grounds 
for theories and approximations. Many well-established methods such as renormalization 
group (RG), high- and low-temperature series expansions etc often failed to deal with these 
systems. Experimentally, data on different frustrated systems show a variety of possibilities: 
first-order or second-order transitions with unknown critical exponents etc. (see reviews in 
Ref. 111). The case of XY and Heisenberg spins on the STA has been intensively studied 
mostly since 1987. There are good recent reviews on the subject (see for example reviews by 
Delamotte et al^). Let us briefly recall here some main historical developments and actual 
situation. In the XY and Heisenberg cases, different materials give different experimental re- 
sults. The anomalous dimension is found negative in many materials and in most numerical 
simulations, the scaling relations are violated and no universality in the exponents was found 
in early simulations. This situation is briefly described in the following. Kawamurai^'^ has 
conjectured by the two-loop RG analysis in d = 3 the existence of a new universality class 
for frustrated magnets. Since then there have been many other perturbative calculations 
with contradictory results.-^^ From 2000, there has been a number of papers by Tessier and 
coworkers —'^ using a nonperturbative RG study of frustrated magnets for any dimension 
between two and four. They recovered all known perturbative results at one loop in two 
and four dimensions as well as for A^ ^ oo. They determined Nc{d) for all d and found 
Nc{d = 3) = 5.1 below which the transition is of first order. However, they found the exis- 
tence of a whole region in the flow diagram in which the flow is slow. As a consequence, for 
A^ = 2, 3, they found pseudo-critical exponents in good agreement with some experimental 
data. This allowed them to account for the nonuniversal scaling observed in XY and Heisen- 
berg frustrated magnets. The only problem in this nonperturbative technique is that the 
Hamiltonian is truncated at the beginning. Giving this fact, we have to be careful about 
its conclusion. As will be seen in this paper, the nonperturbative results are so far well 



confirmed. Early MC results on XY STA have been reviewed by Loison.^° Until 2003, all 
numerical simulations found ambiguous results for this model and never a clear first-order 
transition. A numerical breakthrough has been realized with the results of Itakura^^ who 
used an improved MC renormalization-group scheme to numerically investigate the renor- 
malization group fiow of the Heisenberg and XY STA and its effective Hamiltonian which 
is used in the field-theoretical studies. He found that the XY STA exhibits clear first-order 
behavior and there are no chiral fixed points of renormalization-group fiow for N=2 and 3 
cases. In 2004, Peles et aU^ have used a continuous model to study the XY STA by MC 
simulation. They found evidence of a first-order transition. In 2006, Kanki et al^^, using a 
microcanonical MC method, have found a first-order signature of the XY STA. While these 
recent simulations have demontrated evidence of first-order transition for the XY STA in 
agreement with the nonperturbative RG analysis, all of them suffer one or two uncertain 
aspects: the work of Itakura has used a truncated Hamiltonian, the work of Peles et al has 
used standard MC methods and the work of Kanki et al used a traditional microcanonical 
MC technique. At present, we have a very high-performance technique at hand for weak 
first-order transitions. This is a very good opportunity to test it on the XY STA and to say a 
last word on the nature of the phase transition of this system by using the full Hamiltonian, 
confirming or rejecting the nonperturbative RG and recent MC results. That is the purpose 
of this work. 

We study again here the XY STA with high-resolution MC technique which is very 
efficient specially for weak first-order transition.— Our aim is to try to put an end to the 
controversy which has been lasting for 20 years. We will recall some important numerical 
results in the next section. 

The paper is organized as follows. Section II is devoted to the description of the model 
and technical details of the Wang-Landau (WL) methods as applied in the present paper. 
Section HI shows our results. Concluding remarks are given in section IV. 

II. MONTE CARLO SIMULATION: WANG-LANDAU ALGORITHM 

We consider the stacking of triangular lattices in the z direction. The spins are the 
classical XY model of magnitude S = 1. The Hamiltonian is given by 



H = J5IS.-S, + J'^S,-Sfc (1) 

where Sj is the XY spin at the lattice site i, J2{i,j) indicates the sum over the NN spin pairs 
Sj and Sj in a xy triangular plane, while J2{i,k) indicates that of NN spin pairs between 
adjacent planes. J and J' are in-plane and inter-plane interactions, respectively. We shall 
suppose that J = 1 (antiferromagnetic) and J' = — 1 (ferromagnetic) in the following. 

Recently, Wang and Landau^^ proposed a Monte Carlo algorithm for classical statistical 
models. The algorithm uses a random walk in energy space in order to obtained an accurate 
estimate for the density of states g{E). This method is based on the fact that a flat energy 
histogram H{E) is produced if the probability for the transition to a state of energy E is 
proportional to g{E)~^. 

At the beginning of the simulation, the density of states (DOS) is set equal to one for 
all energies, g{E) = 1. In general, if E and E' are the energies before and after a spin is 
flipped, the transition probability from E to E' is 

p{E^E') = mm[giE)/giE'),l]. (2) 

Each time an energy level E is visited, the DOS is modified by a modification factor 
/ > whether the spin flipped or not, i.e. g{E) -^ g{E)f. In the beginning of the random 
walk the modification factor / can be as large as e^ c^ 2.7182818. A histogram H{E) records 
how often a state of energy E is visited. Each time the energy histogram satisfies a certain 
"flatness" criterion, / is reduced according to / — > v7 and H{E) is reset to zero for all 
energies. The reduction process of the modification factor / is repeated several times until 
a final value /gnai which close enough to one. The histogram is considered as flat if 

H{E) > x%.{H{E)) (3) 

for all energies, where the flatness parameter 0% < x% < 100% controls the accuracy of the 
estimated g{E), with increasing accuracy as x% approaches unity. {H{E)) is the average 
histogram. 

Thermodynamic quantities'*^ can be evaluated using the canonical distribution at any 
temperature T by P{E,T) = g{E)exp{—E/kBT)/Z where Z is the partition function de- 
fined by Z = Y.Eg{E)exp{-E/kBT). 



In this work, we consider a energy range of intereslr^ii^ (i^min, -E'max)- We divide this 
energy range to R subintervals, the minimum energy of each subinterval is -E^in ^"^^ ^ ~ 
1, 2, ..., R, and maximum of the subinterval i is -Emax — -^min + ^A-E', where AE can be chosen 
large enough for a smooth boundary between two subintervals. The Wang-Landau algorithm 
is used to calculate the relative DOS of each subinterval (-Emin; -^max) with the modification 
factor /final = exp(10^^) and flatness criterion x% = 95%. We reject the suggested spin flip 
and do not update g{E) and the energy histogram H{E) of the current energy level E if the 
spin-flip trial would result in an energy outside the energy segment. The DOS of the whole 
range is obtained by joining the DOS of each subinterval (-Emin + ^-^> -^max ~ A^E). 

III. RESULTS 

We used the system size of N x N x N where N = 

12,18,24,30,36,48,60,72,84,90,96,108,120. Periodic boundary conditions are used 
in the three directions. | J| = 1 is taken as unit of energy in the following. 

The energy histograms for two representative sizes A^ = 48 and N = 120 are shown in 
Figs. [Hand [21 respectively. As seen, for A^ = 48, the peak, though very large, does not show 
yet a double-maximum structure. Only from A^ = 90 that the double-peak structure clearly 
appears. This is a sufficient condition, not a necessary condition, for a first-order transition. 
We give here the values of T^ for a few sizes: T^ = 1.458270, 1.457878, 1.457642, 1.457537 
for A^ = 48, 84, 96, 120, respectively. Note that this result is in excellent agreement with 
earlier MC simulationsiiii^ii^ using less sophisticated methods. To explain why standard 
MC methods without histogram monitoring (see for example Ref. |3|) fail to see the first 
order character, let us show in Fig. [3] the energy vs T obtained by averaging over states 
obtained by the WL method for A^ = 120. We see here that even at this big size, the average 
energy does not show a discontinuity as in a strong first-order transition: the averaging over 
all states erases away the bimodal distribution seen in the energy histogram at the transition 
temperature. Therefore care should be taken to avoid such problems due to averaging in 
MC simulations. We note that the distance between to peaks in Fig. [21 i. e. the latent 
heat, is ~ 0.009 in agreement with earlier works.— »^>^ 
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FIG. 1: Energy histograms for A^ = 48 at T^ indicated on the figure. 
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FIG. 2: Energy histograms for N = 84, 96, 120 at T^ indicated on the figure. 

IV. CONCLUDING REMARKS 

We have studied in this paper the phase transition in the XY STA by using the flat 
histogram technique invented by Wang and Landau. The method is very efficient because 
it helps to overcome extremely long transition time between energy valleys in systems with 
a possible first-order phase transition. We found that the transition is clearly of first-order 
confirming therefore recent MC results using less efficient techniques. These results put 




FIG. 3: Energy vs T for iV = 120. 



definitely an end to the 20-year long controversy and lend support to nonperturbative RG 
calculations using an effective average Hamiltonian. 
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